Light diffusion and localization in 3D nonlinear disordered media 
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Using a 3D Finite-Difference Time-Domain parallel code, we report on the linear and nonlinear 
propagation of light pulses in a disordered assembly of scatterers, whose spatial distribution is 
generated by a Molecular Dynamics code; refractive index dispersion is also taken into account. We 
calculate the static and dynamical diffusion constant of light, while considering a pulsed excitation. 
Our results are in quantitative agreement with reported experiments, also furnishing evidence of a 
non-exponential decay of the transmitted pulse in the linear regime and in the presence of localized 
modes. By using an high power excitation, we numerically demonstrate the "modulational instability 
random laser": at high peak input powers energy is transferred to localized states from the input 
pulse, via third-order nonlinearity and optical parametric amplification, and this process is signed 
by a power-dependent non-exponential time-decay of the transmitted pulse. 



Electromagnetic wave localization and diffusion in ran- 
dom media have been recently considered by several au- 
thors. Experimental studies reported on various evi- 
dences of light localization, including the existence of 
a non-exponential long-time tail of the transmitted in- 
tensity [H, 0, H, 0, IE IE 0, S] , and these processes were 
considered in numerous theoretical and numerical papers 
(as e.g. in @, 0, [H, [3, US Q 0, El ) ■ In a purely diffu- 
sive regime, pulsed light impinging on a disordered sam- 
ple is rapidly dispersed while propagating; this implies 
that nonlinear effects are typically negligible. This holds 
true as far as localization processes do not get involved. 
For example, in random lasers (for a recent review see 
localized or extended modes [l8| are excited in a 
disordered medium by means of an optically active ma- 
terial; in this case, the large resonant nonlinear suscep- 
tibility and the fact that long living (high Q) modes are 
excited foster a variety of nonlinear phenomena (as e.g. 



19L 20l 2lj). For non-resonant ultra- fast nonlinearities 
(as x processes of electronic origin [22j), it is in princi- 
ple possible to use high-peak intensity femtosecond laser 
pulses to induce nonlinear effects, even in the presence of 
strong light diffusion. However, the (numerical and theo- 
retical) analysis is enormously complicated by the need to 
include various effects as disorder, nonlinearity, material 
dispersion and to take into account a 3D environment. 
Even in the absence of nonlinearity, no "ab initio" numer- 
ical investigation of 3D Maxwell equations for femtosec- 
ond pulses in disordered media has been reported, to the 
best of our knowledge. Some theoretical investigations 
considered the role of the Kerr effect (i.e. an intensity de- 
pendent refractive index) for a monochromatic excitation 
|23l. |24|. l25l |26| ; however ultra- fast optical nonlinear- 
ity may also provide other mechanisms. We specifically 
consider the gain of non-resonant origin, i.e. the optical 
parametric oscillation (OPO). OPO, in isotropic media 
in the absence of disorder, has been recently considered 



with reference to various devices 27 , 2i| , while in the 
fiber-optics community it is known as the "modulational 
instability laser" [3(|. Here we report on the extensive 
3D+1 numerical analysis of linear diffusion and nonlinear 
amplification processes in a random dispersive nonlinear 
medium. We use a Molecular Dynamics (MD) code for 
generating a three dimensional distribution of spherical 
scatterers and a Finite-Difference Time-Domain (FDTD) 
parallel algorithm for simulating light propagation. The 
MD-FDTD approach not only reproduces the diffusive 
regime in agreement with reported experimental results 
in the linear regime, but also provides clear indications 
of the random OPO process, with analogies to random 
lasers. 

An "ab initio" investigation of the nonlinear light prop- 
agation in 3D disordered media, requires a realistic dis- 
tribution of a colloidal medium and of its nonlinear di- 
electric response. We model the disordered system as an 
assembly of poly-dispersed particles in air. The MD code 
generates the particle configuration; specifically, we con- 
sidered a 50:50 mixture of spheres with diameter ratio 
1.2, interacting by a 100 — 200 Lennard Jones potential 
3l| . The considered sample is composed by 1000 scat- 
terers; once determined the particle positions by the MD 
code, each of them is taken as filled by a dispersive op- 
tically nonlinear medium, which is modeled as described 
below. The particle dimensions are chosen in order to 
simulate typical samples used in light localization and 
diffusion experiments, as e.g. in [l|, |5|. We considered 
diameters for the biggest particles in the mixture rang- 
ing from 200 nm to 400 nm; correspondingly the volume 
filling fraction <f> for our samples, with volume L 3 (L = 4 
/im, of the order of the values considered in [||), varies 
in the interval <fi £ [0.1,0.6]. The optical linear response 
of the particle medium is chosen as the best approxi- 
mation of that of T1O2, used in most of the reported 
experiments, including material dispersion. The latter is 
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modeled by a single-pole Lorentz medium, reproducing 
the data for TiOi reported in [32| . Material absorption is 
also included providing an absorption coefficient Ki = 10 
mm -1 at A = 514 nm. The 3D Maxwell equations, solved 
by the parallel FDTD [33| algorithm, are written as 



V x H = e £ood t E + dtP m 



V x E = -Mo9 t H 

<9 t 2 P + 2 7o 9*P + u; 2 f(P)P = e (e s - £oo ) W 2 E 

The non-resonant nonlinear response is determined by 
the nonlinear Lorentz oscillator (see e.g. [22]). As 
f(P) — 1, eqs.flj) model the linear dispersive medium 
T1O2] by appropriately choosing f(P), the system dis- 
plays a isotropic ultra-fast nonlinearity, as previously de- 
scribed in [27J, including the non- resonant suscepti- 
bility. In the present case, the associated Kerr effect has 
a nonlinear refractive index An = 77.2/, with / the opti- 
cal intensity and ri2 = 10 _18 m 2 W _1 . The parameters in 
TiOi particles are given by (MKS units) : e s — 5.9130, 
£oo = 2.8731 , w = 6.6473 x 10 15 , 7o = 8.9168 x 10 11 
and f(P) = [1 + (P/P ) 2 ] 3/2 with P " 2 = 3.0; out of the 
particles, P = and e s = = 1 (vacuum). 

The code is parallelized using the "MPI" protocol, 
and "UPML" absorbing boundary condition are used in 
the three spatial directions pjjjj . The used discretiza- 
tion is Ax = Ay = Az = 30 nm and At = 0.02 fs. 
The z— propagating input beam is taken with a Gaussian 
TEMoo linearly y— polarized profile, at the input face of 
the disordered medium its waist is wq = 1 /im. The 
input pulse temporal profile is Gaussian, with duration 
to = 100 fs and carrier wavelength A = 520 nm. The 
numerical results have been verified a number of times 
by halving the spatio-temporal discretization. The va- 
lidity of our approach has been confirmed in a variety of 
studies, as e.g. in 27. 34|. 

The first numerical simulations were aimed to inves- 
tigate light diffusion when the nonlinear effects can be 
taken as negligible (low power), for a fixed sample spa- 
tial extension. We considered increasing values for the 
diameter of the spheres in the mixture; in the following 
a is the corresponding value of the biggest particles in 
the considered 50:50 mixture (for the others, the diam- 
eter is a/1.2). This corresponds to increase the average 
density of the medium, as well as the volume packing 
fraction. In order to unveil the onset of diffusion, we 
monitored the time-dependent output intensity, which is 
expected to develop an exponential tail in the diffusive 
regime. This is obtained in two ways: i) we considered 
the electric field E y at one point in the middle of the 
output face of the sample (see figure 1), and we calcu- 
lated the "local intensity" by squaring the field and fil- 
tering with a low pass filter to remove the optical carrier 
(thus mimicking a photodiode); ii) we considered the to- 
tal transmission from the output face of the sample as 
the calculated z— component of the Poynting vector, in- 
tegrated with respect to the transverse x, y coordinates. 



In figure 1 we show the results for the local intensity, 
while increasing the volume filling fraction <f>. In order 
to compare the trailing edge of the transmitted pulses, 
for each signal we rescaled the peak value to the unity 
and shifted the temporal axis to make all the peaks for 
any a coincident. At high packing fraction the onset of 
an exponential trend is clearly visible. Note the residual 
oscillations, that are due to the fact that light is collected 
in a specific point at the output face. A typical far-field 
speckle pattern calculated from our simulation is shown 
in the inset of the bottom panel, and obtained by stor- 
ing the Ey component of the electric field in the output 
section, Fourier transforming with respect to the trans- 
verse coordinates, and averaging over time (the pattern 
is calculated for a continuous wave excitation at A = 520 
nm). 

The intensity modulation in figure 1 disappears when 
considering the overall transmission from the output 
facet, which provides an average information over all the 
speckles and is shown in figure [21 The instantaneous 
flux for the z-component of the Poynting vector is calcu- 
lated and then filtered to remove the reactive component. 
The appearance of an exponential tail in the transmis- 
sion denotes the transition to the diffusive regime. The 
transmitted signal is lit) oc exp[— ir 2 D(t)t/ L 2 ], with an 
instantaneous diffusion coefficient D(t) — > D as t — > 00 
and a sample length L. When the packing fraction in- 
creases, the tail gets longer due to the reduced diffusivity 
D. The latter quantity is shown in Fig. [2{b) vs 4>- The 
results are in agreement with experimentally determined 
values, as reported in literature (e.g. in [3,0,0]). Note 
that, as expected, D decreases as the volume packing 
fraction increases, and at high <f> the diffusivity tends to 
an asymptotic value (for non overlapping particles). The 
inset in Fig. [2jb) shows the calculated instantaneous 
diffusivity for two values of a, which reproduces the re- 
ported trends, including some temporal oscillations [3,0]. 

Next, we consider the diffusive regime for a fixed 
a = 400 nm (<p = 0.6), while increasing the input peak 
power Pi n . In figure [31(a), we show the local intensity for 
two input peak powers; the onset of a non-exponential 
tail is evident, as also shown in the inset reporting the 
intensity auto-correlation function. In the panel (b) we 
show the transmitted intensity (integrated Poynting vec- 
tor) for various input powers. The inset shows the in- 
stantaneous diffusivity that drops at high powers and 
long times, when a threshold value for the input peak 
power is reached. These results foster the interpretation 
of the data in terms of a light-induced localization tran- 
sition, strongly resembling that controlled e.g. by the 
jacking fraction and experimentally observed (as e.g. in 
l[). However, in the present case, the transition is con- 
trolled by the input power, and hence it has a nonlinear 
origin. It will be shown below that it can be explained 
as the excitation of a "modulational instability random 
laser" (or random OPO). 
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Indeed, to unveil the underlying mechanism, we re- 
sorted to the spectral analysis of the electric field in var- 
ious points inside the random medium; figure H] summa- 
rizes the typical result. In panel (a) we show the calcu- 
lated frequency spectra of the electric field E y in the ran- 
dom structure, for different input peak powers. The data 
are displayed in logarithmic scale and vertically shifted 
by an arbitrary amount for the sake of clearness. At 
low power the spectrum of the input pulse is reproduced, 
with the addition of small ripples (due to the vertical 
scale in figure 2^) in correspondence of the existing high- 
Q modes. As the power increases, the bandwidth gets 
larger and regions with depression and enhancement of 
the spectrum appear. At sufficiently high power, a wide 
band excitation is found. This behavior can be explained 
by recalling that in media energy can be transferred 
from various harmonics to others through parametric 
amplification. We stress that the nonlinear effects are 
weakly pertubative, indeed we consider power level just 
above the OPO threshold and the nonlinear index per- 
turbation, due to the Kerr effect, is of the order of 0.001. 
Thus the localization process in figure [3] is explained by 
assuming that energy is nonlinearly transferred to local- 
ized, long-living states. 

In panel (b) we show the time evolution of the spec- 
trum (i.e. the spectrogram calculated by Fourier trans- 
forming windowed parts of the temporal signal) at Pj„ = 
3 kW. In the initial interval t < 1 ps, energy is trans- 
ferred to a wide band around the input pulse spectrum. 
However, at large times only the long-lifetime (high Q- 
factor) modes survive and the spectrum is character- 
ized by several peaks in the region around A = 300 nm 
(v = 0.9 x 10 15 Hz). In order to sustain the fact that 
localized modes put into oscillation, we show in the inset 
of figure 4(c) a snapshot of the electric field E y in the 
sample, taken at P = 2.5 kW and t = 4 ps which clearly 
show the occurrence of localized excitations. 

In addition, we went back to the linear regime to an- 
alyze the localized modes are supported in our virtual 
samples. First, we repeated the pulse transmission ex- 
periment at low power P = 1 nW and central wave- 
length A = 300 nm (a = 400 nm) and Figure 4(c) clearly 
shows an non-exponential tail at t = 2.5 ps, which is not 
present for A = 521 nm (figure 2a) also for t > 3 ps (not 
shown). Then, we recall that light localization is typi- 
cally expected in proximity of depressions of the density 
of states (DOS). For the case of weakly disordered peri- 
odical structures, these are obtained inside the photonic 
band gaps [35| . while, for random dispersions of parti- 
cles, the DOS (which is numerically calculated by shin- 
ing the sample by a ultra-wide band single-cycle pulse 
and spectrally analyzing the transmitted pulse) displays 
some depressions (pseu do-g aps) in proximity of the Mie 
resonances (see e.g. 0, lll|). In our case, the distribu- 
tion of states displays various depressions with the most 
pronounced around 0.9 x 10 15 Hz, as shown in Fig. fjjd). 
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FIG. 1: (Color online) Top panels (a-c): electric field in the 
middle of the output face of the sample at low input peak 
power (Pi„ = 1 nW), for three different particle diameters 
(a — 250,350,400 nm). The insets show the particle dis- 
tribution in the yz (top) and xz (bottom) middle section of 
the sample. Bottom panel (d): local intensity in logarithmic 
scale. The inset shows the calculated speckle pattern. 
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FIG. 2: (Color online) (a) Tail of the transmitted pulse (the 
flux of the z-component of the Poynting vector at the out- 
put face) for various particle diameters. The inset shows the 
transmitted pulse in linear scale, the gray region is the in- 
put pulse, (b) Calculated diffusion constant for various filling 
fractions. The continuous line is a best-fit by an inverse power 
law. The inset shows the dynamical diffusivity for two values 
of the particle diameters. 



The fact that localized modes do exist in proximity of the 
pseudo-gaps, is also confirmed by the plot of the calcu- 
lated Q-factor vs frequency (bold line in Fig. Hid)), real- 
ized by using an harmonic inversion library [36(. Which 
of the localized modes are effectively excited will depend 
not only on the Q-factors, but also on the spatial overlaps 
with the pump modes and with the spatial distribution 
of the nonlinear susceptibility (the g coefficient discussed 
below). The localized modes are a reminiscence of the 
Mie resonances of the isolated spheres, which are strongly 
affected by the presence of adjacent resonators in a disor- 
dered configuration, in analogy to what happen for peri- 
odical mono-dispersed photonic crystal, investigated e.g. 
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FIG. 3: (a) Time dynamics of the local intensity for two dif- 
ferent input peak powers; the inset shows the corresponding 
auto-correlation functions versus the delay r; (b) transmit- 
ted signal for various input peak powers; the inset shows the 
calculated instantaneous diffusivity. 
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FIG. 4: (Color online) (a) Spectral emission for increasing 
peak power; (b) spectrogram displaying the time evolution of 
the spectra for P = 3 kW; (c) transmission of a low power 
pulse at A = 300 nm and a — 400 nm displaying a not- 
exponential tails; the inset shows the field distribution for 
a — 400 nm at high power (P = 2.5 kW) at t — 4 ps in 
the xz middle section of the sample; (d) density of states and 
Q-factor distribution for a = 400 nm. 



in [37(. 

Before concluding we show that the basic features of 
the reported results can be reproduced by a simple the- 
oretical model, based on the coupled mode theory in 
the time domain (CMT) [3^]. Indeed the small-signal 
parametric gain for the four-waves interaction luq + ujq = 
(u) + n) + (Ldo~Q), can be determined by introducing the 
complex amplitudes for the pump mode a (t) = a(u ,t) 
and for the generated modes a±(t) = ci(ujq ± Q,,t). The 
complex amplitudes are normalized in a way such that 
\a(uj,t)\ 2 is the energy stored in the mode at angular 
frequency cj, and we denote with t(cj) the lifetime, and 
Q(ui) = ujt(uj)/2 the corresponding Q-factor. In the un- 
depleted pump approximation the relevant CMT equa- 
tions are written in compact and obvious notation as 



da± 



1 ■ 2 * 

= a± + iga^a^. 



(2) 



tween the mode electric field distributions and the non- 
linear third-order susceptibility. Looking for solutions in 
the form a± = A± exp(af), the real- valued gain a is given 
by 



1 

2^ 



1 

2rl 



(—- — 

V 2t+ 2t_ 



-l 1/2 



\ga. 



l\ 2 



(3) 

It is readily seen that a condition for having amplification 
(a > 0) is given by \ga 2 a \ 2 > l/(r+r_) = lu + u + /AQ+Q-. 
Assuming that for all the involved modes uj = ujq and the 
Q-factor is almost the same, with the exception of those 
localized that display higher values, it is concluded that 
the localized modes (those with not negligible overlap g 
with the pump modes) are put into oscillations at lower 
pump energies (see line at 2.25 kW in figure SJa)). Ad- 
ditionally the asymmetry, at low power, in the spectrum 
(see 0Ji) can be explained by observing from Eq. @ that 
> 1 if Q+ > At higher powers most of the 
modes are put into oscillation and the spectrum mimics 
the DOS. 

In conclusion, the considered disordered, nonlinear and 
dispersive samples support localized modes in proxim- 
ity of Mie resonances. These can be excited by high- 
intensity pulses through optical parametric amplification. 
The overall process can be interpreted as the excitation 
of a "modulational instability random laser" . The main 
difference with standard random lasers is given by the 
fact that an ultra-fast optical process provides the re- 
quired gain and does not need doping of the sample by 
an active media. These results can be generalized to 
other nonlinear amplification processes, like the Raman 
effect, and can find application in all-optical storage of 
light, ultra-fast laser-tissue interaction, soft-matter and 
granular systems. We thank S. Trillo and D. Wiersma 
for fruitful discussions. We acknowledge support from 
the INFM-CINECA initiative for parallel computing. 



The coefficient g is the relevant 3D spatial overlap be- 
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